Implementation of the LDA+U method using the full potential linearized augmented 

plane wave basis 
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We provide a straightforward and efficient procedure to 
combine LDA+U total energy functional with the full po- 
tential linearized augmented plane wave method. A de- 
tailed derivation of the LDA+U Kohn-Sham type equations 
is presented for the augmented plane wave basis set, and a 
simple "second-variation" based procedure for self-consistent 
LDA+U calculations is given. The method is applied to cal- 
culate electronic structure and magnetic properties of NiO 
and Gd. The magnetic moments and band eigenvalues ob- 
tained are in very good quantitative agreement with previous 
full potential LMTO calculations. We point out that LDA+U 
reduces the total d charge on Ni by 0.1 in NiO. 



I. INTRODUCTION 

The limitations of the local density approximation (LDA) 
in calculations of electronic and magnetic ground state prop- 
erties of strongly correlated materials are well knownO. One 
of the simplest methods for going beyond LDA-.is provided 
by so-called "LDA+U" total energy functionaloEl, which has 
proven to give large but reasonable corrections for a number 
of magnetic insulators. (Throughout this paper we under- 
stand LDA to include its generalization to spin polarized sys- 
tems.) The method presumes that an appropriate set of local 
orbitals {<j>m} (such as 3d or 4/) can be identified. Then a 
strong intra-atomic interaction is introduced and treated in 
a Hartree-Fock-like manner; in particular, the intra-atomic 
exchange term is treated without any local density approxi- 
mation. A primary result is the splitting apart of occupied 
and unoccupied states within the {<j> m } shell. 

Technically, the explicit use of local orbitals and an or- 
bitally dependent correction to the energy functional and to 
the LDA effective potential makes it most convenient to im- 
plement the LDA+U within a method using atomic-like or- 
bitals as basis functions. This is true for the linearized muffin- 
tin orbital method (LMTO) in atomic sphere approximation 
(ASA)p or with its full potential versionu, both of which are 
formulated in terms of atom-based orbitals, and most LDA+U 
work has been done with the LMTO- ASA method. Recently 
the implementation of LDA+U in the psaudopotential plane 
wave method was presented by SawadaH by projecting the 
plane waves onto an atomic orbital. 

The LDA+U method has not yet been implemented within 
an all-electron, full potential method whose basis set is not 
based explicitly on local orbitals. The linearized augmented 
plane wave (LAPW) method, which makes no shape approxi- 
mation and is acknowledged to be state-of-the-art in accuracy, 
uses a basis set of plane waves that are matched onto a linear 
combination of all radial solutions (and their energy deriva- 



tive) inside a sphere centered on each atom. Hence the basis 
functions are nothing like individual atomic orbitals. It is im- 
portant to establish that the LDA+U method is not limited 
to certain basis sets, but rather is more widely applicable. 
The full potential aspect may also become important in-del- 
icate situations, such as at the onset of orbital orderingp. In 
this paper we present a simple and numerically efficient pro- 
cedure to combine the LDA+U total energy functional with 
the full .potential linearized augmented plane wave (FLAPW) 
methoda. As an example of its efficacy we apply our imple- 
mentation to study the effects of strong on-site repulsion on 
the magnetic insulator NiO and the ferromagnetic metal Gd 
and compare with previous work on these materials. 

The paper is organized as follows. For the sake of com- 
pletenesaqin Sec. II we recall the basic equations of LDA+U 
methodEla. Then we address the fact that the LAPW method 
is not a local orbital basis set approach. We show that we can 
define the on-site orbital occupation matrix and formulate 
a variational procedure to obtain Kohn-Sham single particle 
equations. We also describe a "second variation" procedure 
that may provide additional insight in analyzing the results. 
In Sec. Ill we present the results of self-consistent total en- 
ergy LDA+U FLAPW calculations for the antiferromagnetic 
insulator NiO and for Gd metal. These results are compared 
with previous work, and additional features are pointed out. 



II. COMPUTATIONAL METHOD 

The essence of the LDA+U method is to identify atomic- 
like orbitals {4>m} (m is the magnetic quantum number, we 
will often suppress other indices) and to treat interactions 
of electrons in these orbitals in a non-LDA manner. The 
LDA+U variational total energy functional takes the for 



E tot (p,n) = E LDA (p) + E ee (h) - E dc (h) 



(1) 



where, E LDA (p) is usual local spin density functional of the 
total electron spin densities p CT (r) (a =t, I)- E ee is an 
electron-electron interaction energy and E dc is a "double- 
counting" term which accounts approximately for an electron- 
electron interaction energy already included in E LDA . Both 
are functions of the local orbital occupation matrix n CT = 
rf mm i of the orbitals {<t> m }- 

A particular form of E ee is taken in accordance to the 
multi-band Hubbard model for d(f) electronsB 
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where the V ee is an effective on-site "Coulomb" interaction 
and < | > indicates an angular integral. 

The electron-electron interaction potential in the atomic 
limit is given by 

< m ly m 3 \V ee \m 2 ,m4 >= a fc (mi, m 3 , m 2 , m 4 )F k (3) 

k 



afc(mi,m 3 ,m2,m 4 ) = 
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where F k are the Slater integrals, \lm > is d(f)-spherical har- 
monic, and n^ ni m2 is the on-site d(f ) occupation matrix in the 
spin-orbital space, which has to be defined with respect to the 
chosen localized orbital basis set. Note that in E ee in Eq.(2) 
the self-interaction in the first term (mi = 7712 = 7713 = m± 
and a = a') is cancelled exactly by the exchange interaction 
in the second term. This is one of the important features of 
the LDA+U approach. The "double-counting" term is taken 
to satisfy an "atomic-like" limit of the LDA total energyu: 

£ dC (n) = ^n(n-l)-^YJn> CT -l) (4) 



»k + c(r) 
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where ui , 1I1 are the radial wavefunction and its energy deriva- 
tive at an appropriate energy Ei for angular quantum number 
/ (recall < ui\ui >= 1, < ui\ili >= 0), k is a k-point in BZ, 
and G is any reciprocal lattice vector. 

The solution of the Kohn-Sham equations is a set of eigen- 
functions 



(8) 



where b stands for band index and the sum runs over recip- 
rocal lattice vectors. The charge density is 



(9) 



where /k,i, is the occupation of the state. In the i-th sphere, 
from Eqs. (8) and (9) the spin density is 



P"( T ) = Ek,b /kf>E G ,G' C k+G' C k+G >< 
E K+G' U V 00 + 6 k+G< 00] 



(10) 



where the Hubbard U and exchange J constants are given by 
1 1 



U = 
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J = U 



21(21 + 1) ^ 



[< mi, mz\V ee \mi, m 3 > 



- < mi,m 3 |V e<: |m 3 ,mi >] . 

In Eq.(4), n" = Tr n a , and n = ^ ct n CT is a total d(f) on- 
site occupation numbers. The energy Eq.(l) is a functional 
with respect to the spin density 

P°(r)=J2fi*r(rW(r) 0) 

i 

and the local orbital occupation matrices n" . 

A. Density matrix definition within the LAPW basis 

It is rather straightforward to choose an orbital basis for 
Eq.(3) in the case of the LM.TO method as an orthogonal, 
atomic-type LMTO basis settl In the case of the LAPWB 
method there are no naturally defined atomic-like basis or- 
bitals. Therefore, we will use a projection of augmented plane 
waves onto radial functions inside muffin-tin spheres to define 
the density matrix for Eq.(3). 

Let us recall some basic features of LAPW method (we 
will use the notations of Singho). The LAPW basis func- 



tion 0k+G( r ) i s a planewave exp ^i(k + G) ■ rj outside MT 

sphere, and inside sphere i (vi = r — R;) it is a linear combi- 
nation of all radial functions 



[flk+euf (n) +&L+ G wf00] Yim(n)Y; m ,(n). 

Taking the elements of the density matrix from the projection 
of the wavefunctions onto the Yi m subspace gives 
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Such a projection gives the desired result (i.e. in agreement 
with other implementations) in the atomic limit, and its di- 
agonal elements give what are commonly thought of as local 
orbital occupation factors. 

This expression indicates that the density matrix is not 
necessarily diagonal in I but this detail is not expected to be 
important. Keeping only the I' — I part and suppressing the 
I index on the density matrix, after radial integration Eq.(ll) 
becomes 
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Finally, one can rewrite Eq. (12) as 
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B. Variational principle. 

Using Eq.(13) and orthonormaAjty of the $'s one can min- 
imize Eq.(l) with respect to $* ,<T Lj: 



5*. 



= 



(14) 



It gives the set of Kohn-Sham equations, which can be written 
in the form 



(-V 2 + ^Vt(r))<&r(r) + 

E ^ 



(15) 
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where an effective potential acting on the Yj m subspace (the 
d(f)-states) is 



(16) 
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Eqns (15)-(17) indicate specifically how the additional poten- 
tial that arises in the LDA+U method acts on the projections 
of the wavefunctions <E> onto the Yi m subspace. 



C. Use of "second variation" procedure. 

The Eq.(15) can be solved directly using the LAPW basis 
set Eq.(7). For large systems it can be solved more efficiently 
using a "second variation" procedure. We define an auxiliary 
orthogonal basis set ^' CT (r) as the solutions of LDA band 
Hamiltonian 



(-V 2 + ViV 4 (r))*r (r) = e 6 ' CT *r(r) 



(18) 



The functions are somewhat like the LDA eigenfunctions, 
but they are not identical because it is not the LDA den- 
sity that goes into V£ DA . Expanding the functions $ via >]/ 
(suppressing k and shortening notation) 



(19) 



Eq.(15) is then transformed to the following, 

J2 Eidj l* 3 >+E^ E V ™™' [< U l Y lm'\* j > X (20) 
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and using < "J;!^ >= Sij we obtain the "second variation" 
Hamiltonian Hjjid^, = e^d* ) of the form: 
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Now, given some plane wave cutoff as usual, the self- 
consistent solutiep of Eq.(15) is performed as follows, (i) 
Solving Eq.(18)E-3 for the given LDA effective potential we 
obtain the entire orthogonal basis set {"I"}; (ii) within this 
basis set we solve Eq.(21) to obtain the coefficients dj of the 
LDA+U wavefunctions in terms of the auxiliary wavefunc- 
tions in Eq.(19); (iii) the new LDA+U wavefunction is then 
projected back to LAPW basis set (cf., Eq.(7)) and the usual 
procedure to calculate charge and spin densities Eq.(10)H-J 
and occupation matrices Eq.(12) is then used to achieve self- 
consistent solution. 



D. Total energy evaluation 

After the self-consistent solution of Eq. (15) the LDA+U 
total energy can be calculated from Eq. (1), using: 



E L 



(p)=T s (p) + / drv ext (v)p(r) + E H {p) + E xc (p) 

(22) 



where T s is the kinetic energy for non-interacting electrons 
with density p, v ext is an external potential (including the 
interaction with nuclei), Eh is the Hartree energy and E xc 
is an exphange-correlation energy. It is convenient to use, 
as usualE-3, the eigenvalue sum, which includes T s and the 
interaction with nuclei exactly but miscounts other terms. 
A so-called "double-counting" correction Ef^ DA leads to the 
correct energy. 

For our LDA+U calculations, the eigenvalue sum contains 
additional terms (involving U and J) whose contributions to 
the energy are miscounted. When these terms are corrected, 
the expression for the energy is: 



E tot (p,h) = E £i ~ e lda(p) ~ E d L c DA+u {h) 



(23) 



where, Ef c DA is a "double-counting" LDA correction, and 
Elda+u i s a "double-counting" correction to the eigenvalue 
sum in LDA+U calculations. This correction is expressed 
in terms of electron-electron interaction energy Eq. (2) and 
"double-counting" term Eq.(4) as 

Et c DA+u (p, n) = E^{n) - E dc (n) - 1(17 - J)n (24) 

where, the last term appears due to a specific form of the 
LDA+U "double counting" energy Eq. (4) with an explicit 
removal of self-interaction. 
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The LDA+U expression for the energy (Eq. (1)) has been 
used very sparingly except to generate the Kohn-Sham equa- 
tions which lead to the self-consistent density and orbital oc- 
cupation matrices. Partly this has been due to the imple- 
mentation of the LDA+U in the LMTO-ASA method, which 
makes severe approximations to the shape of the potential and 
density and for the variational freedom of the wavefunctions, 
which translates into questionable ability to resolve energy 
differences. Another feature of the LDA+U energy functional 
that has not been emphasized is that it explores a regime of 
Elda[p] which is not understood, by evaluating Eld a for 
densities that do not minimize it. This procedure deserves 
scrutiny, so we will give attention to the energy differences in 
the following section. 

III. RESULTS 

As representative systems to illustrate the numerical pro- 
cedure we have proposed, we choose the antiferromagnetic in- 
sulator NiO and Gd metal, which has a ferromagnetic ground 
state. For both of them LDA fails to provide the correct elec- 
tronic and magnetic ground states due to significant electron 
correlation effects. 



primarily Ni e g states in the conduction band as opposed 
to the LDA calculations where a small band gap (0.4 eV) is 
formed between mainly Ni d states. There is a large downward 
shift of majority Ni e 9 states (cf. Fig. 2) and the resulting 
energy difference between e g and e 9 states («i-40.8 eV) is 
not far from model GW calculations (w 9 eV)lij. For the 
symmetry of the Ni ion, the occupation matrix is specified by 
two numbers (for each spin) , the e 9 and t 2g occupancies. The 
LDA+U correction results in a decrease in the Ni d occupation 
by 0.12 electrons, with an increase in spin majority charge of 
0.2 electrons (to essentially 100%) more than compensated by 
a decrease by 0.3 electrons in spin minority charge (nearly all 
of it e g ). This difference of d charge is potentially important 
for the understanding of correlation in magnetic insulators. 

Since we take into account only the muffin-tin part of 
LDA+U correction, it is important to analyze numerically 
the stability of our results with respect to the choice of the 
MT radii. We have varied the sphere radii for both Ni and O 
atoms in the unit cell and found that these changes have little 
effect on the calculated values of magnetic moment and band 
gap (cf. Table |n|)), due to the fact that the sphere contains 
all but a few percent of the Ni d states. 



B. Gd 



A. NiO 

lda-asaIS calculations of antiferromagnetic (AFII) NiO 
gives the value of the band gap and magnetic moment in dis- 
agreement with experimental data (cf., Table [i]). While the 
band gap requires an energy dependent self-energy to obtain 
rigorously and therefore is not a true test of LDA+U correc- 
tions, the magnetic moment is a fundamental test. Generally, 
it is expected that an improved energy functional will also 
give a band gap in closer agreement with experiment than 
the very poor LDA value. Our FLAPW LDA calculations 
give somewhat larger values of both moment and gap than 
does the ASA calculation (Table |j|), reflecting the approxi- 
mations in the shape of the density and potential, and the 
less flexible basis functions in the ASA. For NiO, we choose 
the lattice constant as 7.927 a.u. in accordance to the exper- 
imental datatj. For self-consistency, 110 special k-pointst£l in 
the irreducible part of the AFII BZ were used, with Gaussian 
smearing of 2 mRy of the Fermi energy to promote conver- 
gence. The sphere radii Rmt — 2-0 a.u. and Rmt = 1-8 a.u. 
were used, and Rmt x K ma x = 6.6 determined the basis set 
size. Literature values U = 8.0 eV and J = 0.95 eVul were 
used to obtain the values of Slater integrals in Eq.(3) (Fo = 
8.00 eV, F 2 = 8.19 eV, F 4 = 5.11 eV). 

The calculated values of spin magnetic moment and band 
gap are shown in Table |, are within about 8% of recent full- 
potential LMTO based LDA+U calculations!! The calculated 
projected densities of states (PDOS) for Ni d states and O p 
statep are shown in Fig. 1 and are very similar to these of 
Reffl. The t 2g and e g PDOS for Ni d states for LDA+U 
are compared with those for LDA in Fig. 2. The LDA+U 
calculations give a charge transfer type band gap (~ 3.4 eV) 
between mainly O p states at the top of valence band and 



A detailed analysis of ground state and magnetic prcpi 
erties of ferromagnetic Gd has been presented by SinghJii 
It was shown that Gd 4f-electrons have to be considered as 
"band" electrons instead of "core" to obtain at least qualita- 
tive agreement with experimentally observed Fermi-surface. 
(The "core" treatment consists of treating the 4f states as 
corelike, nondispersive bands that do not hybridize with the 
valence bands, and removing the 4f character from the valence 
electron basis set.) Heiniemann and Temmerman found with 
the LMTO-ASA methodta that using GGA rather than LDA 
led to a ferromagnetic state being favored over the AF state. 
The more recent full potential LMTO calculations of Harmon 
et alx-i did not confirm this conclusion. Instead, it was shown 
that only LDA+U gives a correct ferromagnetic ground state 
for Gd. 

For Gd, we chose the lattice constant as 6.8662 a.u. and 
c/a ratio as 1.58.7 in accordance with the experimental crystal 
structure datallj. Here, 150 special k-points in irreducible part 
of the hexagonal BZ were used, again with Gaussian smearing. 
The values Rmt = 3.2 aaj~ and Rmt x K max = 9.6 were used. 
Literature values of RefO U =6.7 eV and J — 0.7eV were 
used to obtain the values of Slater integrals in Eq.(3) (Fq = 
6.70 eV, F 2 = 8.34 eV, F 4 = 5.57 eV, F 6 = 4.13 eV ). 

The spin magnetic moment for FM and AFM Gd are shown 
in Table llll] for both LDA and LDA+U calculations. The 



LDA result was in nearly perfect agreement with experiment, 
and the 2% increase given by LDA+U does not degrade the 
agreement much. The energy difference E(AFM) — E{FM) 
is negative in LDA calculations and positive in LDA+U cal- 
culations. Our raeults therefore confirm quantitatively the 
conclusion of Ref.til that the correct FM result is only pre- 
dicted by LDA+U. The LDA and LDA+U partial f-densities 
of states are shown in Fig. 3. There is a large (4.5 eV) 
shift of spin-majority f-states below the Fermi level and re- 
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moval of the unoccupied spin-minority f-states from close 
vicinity of the Fermi level due to an upward shift of 1.5 eV. 
The enhancement of an exchange splitting can be understood 
as an increase of an effective exchange splitting parameter 
I — (U + 6J)/7 = 1.5 eV in mean- field Hubbard model in 
comparison with LDA (I — J — 0.7 eV) due to Hubbard 
U. The exchange splitting (IM) for f-states is then increases 
from 5 eV with LDA to 11 eV with LDA+U (cf. Fig.3). A 
more complete discussion of the LDA+U description of Gd 
will be given elsewhere. 



IV. CONCLUSION 

To summarize, we have presented a straightforward and 
numerically efficient procedure to combine the LDA+U to- 
tal energy functional with the FLAPW method. This im- 
plementation is all-electron and includes no shape approxi- 
mations. In spite of a modest limitation - we use the only 
"muffin-tin" part of the LDA+U correction - the method 
works well for representative systems NiO and Gd for the 
electronic spectrum and predicts the correct magnetic or- 
der. It alsp-pan be easily extended to incorporate spin-orbit 
interactiontj. Moreover, a similar numerical procedure can 
be used to implement the "dynamic mean field" parametriza- 
tion within the-ffLAPW method as a more realistic model for 
the self-energyt-3 of metallic strongly correlated systems. 

We have emphasized the eigenvalue shifts in NiO and in Gd 
that result from the strong repulsion in the LDA+U method. 
These shifts are probably important in obtaining an improved 
description. Hovewer, it should be kept in mind that the 
LDA+U method is based on the energy functional for the 
electronic ground state, and that eigenvalues are not the most 
fundamental quantaties for judging the success of the LDA+U 
method. It is satisfying that the LDA+U energies for Gd lead 
to the observed magnetic ground state. 

We are grateful to I. V. Solovyev and A. J. Freeman for 
helpful comments, and R. Weht and CO. Rodriguez for dis- 
cussion. This research was supported by National Science 
Foundation Grant DMR-9802076. 
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TABLE I. The spin magnetic moments (M s in /xs) and 
band gap (E g in eV) for NiO. The increases in both quantaties 
are about the same in the FLAPW and LMTO-ASA methods, 
althrough the individual value differ due to approximations 
made in the LMTO-ASA method. 





FLAPW 


lmto-asaU 


FLMTol 


Exp.U 




LDA LDA+U 


LDA LDA+U 


LDA+U 






1.186 1.687 


1.0 1.56 


1.74 


1.64-1.70 


E 3 


0.41 3.38 


0.2 3.1 


3.4 


4.0 - 4.3 



TABLE II. The spin magnetic moments (M s in /is) and 
band gap (E g in eV) for different "muffin-tin" radii (a.u.) for 
Ni and O atoms in NiO. The small differences are entirely 
unimportant. 



Ni 


2.0 


2.1 


2.1 


O 


1.8 


1.8 


1.7 


M s 


1.687 


1.700 


1.700 


E g 


3.38 


3.38 


3.38 



TABLE III. The spin moments (in /x_b) and the total en- 
ergy difference (in meV) between antiferromagnetic (AFM) 
and ferromagnetic (FM) phases for the bulk Gd. 





LDA 


LDA+U 


Exp. 




FM 
AFM 


7.659 
7.246 


7.818 
7.437 n 


7.63 




E(AFM-FM) 
per atom, meV 


LDA+U 
63 


LDA+Uiil 

68 


LDA 
-2 


LDA-ASA0 
GGA-ASA0 



+ 1.4 



4 - 
3 - 
2 - 




_g I i I i I i I i I i I i I i I i I i I i I i I i I i I i 
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FIG. 1. The partial DOS for Ni d-states (filled) and for 
O p-states for NiO. Majority is plotted upwards; minority is 
plotted downward. 
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FIG. 2. The partial DOS for Ni t 2g (filled) and e g states 
for LDA+U (a) and LDA (b) for NiO. The increase in the 
exchange splitting of the Ni d states due to LDA+U is dra- 
matic. 
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FIG. 3. The total DOS (full line), the partial DOS for Gd 
f-states for LDA+U (filled) and LDA (dotted) for the ferro- 
magnetic Gd. The 4f exchange splitting is increased from 5 eV 
to 11 eV by LDA+U. Majority is plotted upwards; minority 
is plotted downward. 
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